Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.
N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
mkdir -p ~/EvaluationM4M5/FASTQ
mkdir -p ~/EvaluationM4M5/CLEANING
mkdir -p ~/EvaluationM4M5/MAPPING
mkdir -p ~/EvaluationM4M5/QC
cd ~/EvaluationM4M5
tree ~/EvaluationM4M5
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
module load sra-tools
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir FASTQ
Combien de reads sont présents dans les fichiers R1 et R2 ?
cd FASTQ
cat SRR10390685_1.fastq | echo $((`wc -l`/4))
cat SRR10390685_2.fastq | echo $((`wc -l`/4))
Les fichiers FASTQ contiennent 7066055(R1) + 7066055(R2) = 14132110 reads.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
cd ../MAPPING/
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
Quelle est la taille de ce génome ?
zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | awk '/^>/{if (l!="") print l; print; l=0; next}{l+=length($0)}END{print l}' |paste - -
La taille de ce génome est de 4215606 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
Combien de gènes sont connus pour ce génome ?
zgrep -c $'\tgene\t' GCF_000009045.1_ASM904v1_genomic.gff.gz
4448 gènes sont recensés dans le fichier d’annotation.
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
cd ../
module load fastqc
gzip FASTQ/SRR10390685_1.fastq
gzip FASTQ/SRR10390685_2.fastq
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8
La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
car R2 comme le montre a plus de base dans la zone rouge que R1
Lien vers le rapport MulitQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
car il était nécessaire d’enlever “adapter content”
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
La profondeur de séquençage est de : 2.119142 G (total bases before filtering fastp) / 4215606 (ref genome size) = 502.7 (rather good)
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
module load fastp
srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json CLEANING/fastp.json
Les paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| mean quality | >= 30 | to avoid low quality reads |
| sliding window | 8 | sliding window from 3’ extremity to 5’ extremity |
| length of the read | >= 100 | to avoid short reads |
Ces paramètres ont permis de conserver 13.554096M reads pairés, soit une perte de 95,90992427882319% des reads bruts.
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
cd MAPPING/
GCF_000009045.1_ASM904v1_genomic.fna.gz > GCF_000009045.1_ASM904v1_genomic.fasta
module load bwa
srun bwa index GCF_000009045.1_ASM904v1_genomic.fasta
srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fasta ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685_on_GCF_000009045.1_ASM904v1.sam
module load samtools
srun --cpus-per-task=8 samtools view --threads 8 SRR10390685_on_GCF_000009045.1_ASM904v1.sam -b > SRR10390685_on_GCF_000009045.1_ASM904v1.bam
srun samtools sort SRR10390685_on_GCF_000009045.1_ASM904v1.bam -o SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
srun samtools index SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
srun samtools idxstats SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.idxstats
srun samtools flagstat SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.flagstat
samtools faidx GCF_000009045.1_ASM904v1_genomic.fasta
samtools faidx CP031214.1.fasta
Combien de reads ne sont pas mappés ?
13.57 M - 12.83M = 0.74 M reads ne sont pas mappés. Data from MulitQC report
13571369 - 12826829 = 744540 reads ne sont pas mappés. Data from flagstat and idxstats
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
Pour faire cette tâche, j’ai extrait des informations sur trmNF du fichier gff, puis fait manuellement trmNF.lit le fichier. Les fichiers gff et bed utilisent différents systèmes de coordonnées donc, les coordonnées de bed pour le gène = coordonnées gff - 1
zgrep trmNF MAPPING/GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '$3=="gene"' > trmNF.gff3
module load bedtools
bedtools bamtobed -i MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.bed
bedtools intersect -a MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.bed -b trmNF.bed -f 0.5 > result_trmNF_bed.txt
grep -c 'NC_000964.3' result_trmNF_bed.txt
2797 reads chevauchent le gène d’intérêt.
Generation of MulitQC report
export LC_ALL=en_US.UTF-8
export LANG=en_US.UTF-8
module load multiqc
srun multiqc -d . -o .
Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our three affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France Institute Curie, Paris, France